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Using a modified spin-wave theory which artificially restores zero sublattice magnetization on 
finite lattices, we investigate the entanglement properties of the Neel ordered Ji — J2 Heisenberg 
antiferromagnet on the square lattice. Different kinds of subsystem geometries are studied, either 
corner-free (line, strip) or with sharp corners (square). Contributions from the no = 2 Nambn- 
Goldstone modes give additive logarithmic corrections with a prefactor nG/2 independent of the 
Renyi index. On the other hand, corners lead to additional (negative) logarithmic corrections with 
a prefactor 1^ which does depend on both no and the Renyi index q, in good agreement with scalar 
field theory predictions. By varying the second neighbor coupling J2 we also explore universality 
across the Neel ordered side of the phase diagram of the Ji — J2 antiferromagnet, from the frustrated 
side 0 < J2/J1 < 1/2 where the area law term is maximal, to the strongly ferromagnetic regime 
— 1 with a purely logarithmic growth Sq = ^InN, thus recovering the mean-field limit 

for a subsystem of N sites. Finally, a universal subleading constant term 7°'^'^ is extracted in the 
case of strip subsystems, and a direct relation is found (in the large-S limit) with the same constant 
extracted from free lattice systems. The singnlar limit of vanishing aspect ratios is also explored, 
where we identify for 7°'^^ a regular part and a singular component, explaining the discrepancy of 
the linear scaling term for fixed width vs. fixed aspect ratio subsystems. 

PACS numbers: 02.70.Ss,03.67.Mn,75.10.Jm,05.10.Ln 


I. INTRODUCTION 

Entanglement properties of interacting quantum spin 
systems have recently attracted a lot of interest. In par¬ 
ticular, great attention is paid to the universal informa¬ 
tion carried by bipartite entanglement measures such as 
the Renyi entanglement entropies (EEs) defined by 

S', = Y^lnTr(/3o)'', (1.1) 

where pn is the reduced density matrix of a given subsys¬ 
tem fl (see Eig. 1) computed in the ground-state wave- 
function. Note that the special limit of g —)■ 1 cor¬ 
responds to the standard von Neumann EE given by 
Si = —Trpfi In pn and is always implicitly understood 
whenever we refer to g = 1. As a general result, at T = 0 
the Renyi EEs follow an area law^’^ in dimension d 

Sq = aqL^-^ + ■■■ (1.2) 

where L‘^~^ is the size of the boundary between sub¬ 
system n and the rest, and the ellipses are subleading 
corrections. Such corrections have been shown to carry 
universal information about topological order^^®, or the 
presence of Nambu-Goldstone modes associated to the 
breaking of a continuous symmetryIn the latter 
case, Metlitski and Grover (MG)^ have derived the fol¬ 
lowing analytical expression in the case of smooth bound¬ 
aries (no corner), as for instance depicted for d = 2 in 
Fig. 1 (a) for L X ^ strip subsytems: 



where ps is the stiffness, v the velocity of the tiq Nambu- 
Goldstone modes, and 7 °“''^ a universal geometric con¬ 
stant. In the case of subsystems having sharp corners, as 
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Figure 1. Schematic picture for the L x L square lattice Ji — 
J2 antiferromagnet on a torus with 2 types of entanglement 
bipartitions: (a) rectangular strip of extent L x I with no 
corners and (b) square of extent i x £ with four 7r/2 corners. 
In all this work, we assume periodic boundary conditions in 
both directions. 

depicted in Fig. 1 (b), it is expected that^: 

Sq = aqL^-^ + f\n[^L^-^) 

+ +7”^ (1.4) 

where a is a non-universal length scale, and the corner 
contribution depends on no, the Renyi parameter q, and 
the number of corners c of angle pc- The contributions 
IqiPc) from each corner come from the (free) Goldstone 
modes and can be computed, following the work of Casini 
and Huerta^ ^ on scalar field theory, by the numerical 
solution of a set of non-linear differential equations, valid 
for Pc e [0, tt ] {lq{p) = lq{2TT - p)) and g S N \ {1}. 

Previous works have explored the scaling of the entan¬ 
glement entropy in ground-states of systems that break 
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continuous symmetries in the thermodynamic limit. Sub¬ 
leading logarithmic corrections arising from the Gold- 
stone modes have been observed in quantum Monte Carlo 
simulations of finite spin systems^^^^^, even though the 
prefactor of this correction did not perfectly agree with 
the prediction nG/ 2 , until a very recent large-scale, low- 
temperature quantum Monte Carlo (QMC) investigation 
by Kulchytskyy et al.^ for the 2d XY model and q = 2. 
Logarithmic corrections have also been observed in finite- 
size SW calculations® (similar to the ones presented in 
this manuscript), but not with a high-enough precision to 
again ascertain the prediction, except for the case a line¬ 
shaped subsystem for which the prefactor nc/2 could be 
recovered assuming further subleading corrections^^ (see 
also the recent work Ref. 15). The existence of logarith¬ 
mic corrections have also been discussed based on a phe¬ 
nomenological picture of the tower of low-lying states in 
the symmetry-broken phase of antiferromagnets^^. Log¬ 
arithmic corrections due to corner contributions have 
on the other hand been identified and calculated pre¬ 
cisely in free lattice systems^^, broken continuous sym¬ 
metries systems^^ as well as for various critical points 
using QMC, cluster expansions or tree tensor network 
techniques^^^^'^d6-24^ ^ recent work^®, predictions for 
the universality of corner contributions in various theo¬ 
ries are also provided. Finally, Kulchytskyy et al.^ could 
also compute with QMC the subleading constant correc¬ 
tion 7 ”'^ in the 2d XY model, finding a good agreement 
with the prediction of MG in Ref. 7. 

In this paper, we provide a systematic high-precision 
study of the universal nature of three subleading terms of 
the Renyi EE appearing in Eq. (1.4) for a generic model of 
quantum antiferromagnetism in two dimensions (d = 2 ). 
This is achieved using a large-S” semi-classical approach, 
the modified linear spin-wave (SW) theory, where the 
rotational SU(2) symmetry, while practically broken, is 
artficially restored for finite size systems^®’^^. We focus 
on the Ji — J 2 spin-S' antiferromagnet defined on a bi¬ 
partite Lx L square lattice by the following Hamiltonian 

H ^ S, • S, + J 2 ^ s, • s, 

{ij) {{ij)) 

+ (1.5) 

i 

where S are spin-5' operators, interactions act between 
nearest neighbours (ij) and second nearest neighbours 
{{ij)) along the diagonals of a square lattice (see Fig. 1), 
and h is an external staggered field. We impose peri¬ 
odic boundary conditions in all directions. At h = 0 this 
model spontaneously breaks the SU(2) symmetry at zero 
temperature in the thermodynamic limit, and displays 
Neel order for J 2 < J|, with Jtj —t Ji/2 for 5 —>■ 00 ^®. 
The restoration of zero sublattice magnetization in finite 
systems is made possible by tuning the small staggered 
field h*{L) such that on any site (5f) = 0. As first done 
in Refs. 8 and 10, this allows to correctly compute Renyi 
EEs on finite systems. Here we make a systematic and ex¬ 
tensive study across the full Neel regime —00 < J 2 < 


for various subsystem shapes and sizes in order to char¬ 
acterize contributions form (i) Nambu-Goldstone modes, 
(ii) corners, (iii) frustration effects J 2 /J 1 > 0 , and (iv) 
geometric effects appearing through the universal con¬ 
stant 7 ”® in Eq. (1.3). 

Let us briefly summarize our main results. Us¬ 
ing a large-5 approach, we have numerically extracted 
the three subleading corrections in the scaling of EEs 
Eq. (1.4) with no = 2 for SU(2) antiferromagnets. Uni¬ 
versality has been tested in the entire Neel ordered regime 
of the Ji — J 2 Heisenberg model Eq. (1.5) for various 
5, even in the frustrated regime where QMC is inap¬ 
plicable. In the case of subsystems having sharp corners, 
small negative corner terms are found, in perfect agree¬ 
ment with the predictions by Casini and Huerta for free 
scalar fields^The non-universal area-law term has also 
been studied as a function of the second neighbour cou¬ 
pling, showing remarkable behaviors both in the mean- 
field limit (—J 2 /J 1 S> 1 ) where it vanishes, and close 
to the frustrated critical point where the area law 
prefactor strongly increases, while log corrections due to 
Nambu-Goldstone modes are still present. Furthermore, 
the additional geometric constant y”®, which depends on 
the subsystem aspect ratio £/L, is extracted for various 
Renyi indices, and a simple relation with the free scalar 
field result is derived. We have also explored the limit of 
vanishing aspect ratios where a non-trivial slow singular 
behavior shows up as 7 °''‘^(£/L <C 1 ) —t — 00 . 

The rest of the paper is organized as follows. In Sec¬ 
tion H we start by recalling the modified SW formalism 
for the Ji — J 2 spin-5 antiferromagnet, and how it can 
be used to compute the Renyi EEs. We then turn to the 
results for EEs in Section HI where we discuss several as¬ 
pects: we first describe numerical diagonalization results, 
which can be conveniently performed up to subsystems 
of < 10 ® sites, for various shapes of subsystems including 
strips (Sec. HI A and Eig. la) and squares (Fig. lb), with 
a particular focus on the corner contributions (Sec. HI B) 
and their dependence on the Renyi parameter q. In Sec¬ 
tion HI C the dependence on the second neighbour cou¬ 
pling J 2 is studied, focussing on the non-universal area 
law prefactor a^. In Section IV we discuss the constant 
term y”® which is compared to the field-theory predic¬ 
tion of MG in Sec. IV A. An interesting connection to 
the free scalar field result is achieved in Sec IV B. We 
further explore the singular limit of vanishing aspect ra¬ 
tios in Sec. IV C using quasi-analytical results for single 
and double line subsystems where translation symmetry 
inside the subsystem allows to get an explicit expression 
for Sq. Finally we summarize and discuss our results in 
Section V. Details of spin-wave calculations are provided 
in Appendix A, analytical results for the mean-field limit 
— J 2 IJ 2 1 are presented in Appendix B, and an ana¬ 
lytical derivation for one-dimensional subsystems is given 
in Appendix C. 


3 


II. MODIFIED SPIN-WAVE APPROACH 


A. Dyson-Maleev transformation and Bogoliubov 
diagonalization 


We use the Dyson-Maleev formalism^®’^° to map spin 
operators onto bosonic ones. For sites on sublattice A of 
the square lattice: 

s: = s- b\b^ , 5+ = (25 - b%)b^ , S- = bl (2.1) 
and for the sublattice B: 

Sf =b\h^- S, 5+ = -bl{2S -blb^), 5" = -b^. (2.2) 

Truncating at 1/5 order, the Ji — J 2 Hamiltonian 
Eq. (1.5) becomes (up to a constant) 

^ = + (2.3) 

i 

+ Y, SJ2 {bib^ + b]b,) - Y sji {b,b^+b]bi). 

(ij) (iJ) 

After a Fourier transformation, it reads 

n = ^ Ak(5|,6^ + 6L^6_k) + Hk(6|,6Lk + M-k). (2-4) 

k 

with 

Ak = 25 J 2 cos kx cos ky 25( Ji — J 2 ) + ^ (2.5) 

Hk = —5 Ji [cos kx + cos ky] . (2.6) 

The quadratic part of the above Hamiltonian can be di¬ 
agonalized via a standard Bogoliubov transformation: 

6k = MkOk - Wka^k bl = Mkttj^ - WkCt-k- (2.7) 


The quasiparticle operators and satisfy bosonic 
commutation relations provided = 1, and diago¬ 

nalize (2.4) if 



( 2 . 8 ) 

(2.9) 


In terms of Bogoliubov quasi-particles, the J 1 — J 2 Hamil¬ 
tonian takes the simpler form 


% = Yh ^kakCCk + constant, (2-10) 

k 


with the SW excitation spectrum Hk = 2 -\/A^(--/bJ (this 
spectrum is illustrated in Appendix A). In the vicinity of 
the two minima at k = and (0,0), the dispersion 

is linear, with a velocity 


'^sw — 


which is defined only on the AF side J 2 < Ji/2. The SW 
spectrum and velocity are illustrated in Fig. 13 and 14 
of Appendix A. 

In the thermodynamic limit, the continuous SU(2) 
symmetry of the original Ji — J 2 Hamiltonian can be 
spontaneously broken, with the two associated Nambu- 
Goldstone modes at k = (tt, tt) and (0,0). The cor¬ 
responding staggered magnetization order parameter is 
given at the 1/5 order by 


TOAF= lim lim I (5/) I 

h—>-0 N—>-oo 



1 

Stt^ 



( 2 . 12 ) 


In Appendix A, this expression is evaluated numerically 
to obtain the range of parameter space where Neel order 
is expected from this SW treatment. 


B. Spin-Wave theory for finite size systems 


The above SW approach assumes a classical ordered 
state as a starting point. This does not allow for a cor¬ 
rect study of finite size effects since the spin rotational 
symmetry has to remain unbroken on finite-size lattices. 
In order to repair this, adding a staggered magnetic field 
to the quantum antiferromagnet allows to artificially re¬ 
store zero sub-lattice (SW-corrected) magnetization, as 
originally proposed in Refs. 26 and 27. This will turn 
crucial to capture the subleading scaling terms in the en¬ 
tanglement entropy. 

In this approach, one imposes that for any given finite 
size sample (5f) = 0 Vf, which yields a staggered field 
h* such that the number of Holstein-Primakoff bosons 
(n) = 5. This leads to 


Ak(h*) 

V^2kM 


= N{2S+1). 


(2.13) 


This regularizing field is very small and scales rapidly 
to zero with the system size®. Indeed, one can rewrite 
Eq. (2.13) as follows 


iV(25 + l)- ^ 

k^^ko 


Ak(/i*) 

Hk(/i*) 


o ^ko(fe*) 

^2ko(/i*)’ 


(2.14) 


where ko = (0, 0) and (tt, tt) are the singular modes where 
the dispersion vanishes in the absence of staggered field. 
The contributions from these two modes, divergent in the 
limit /i* —>■ 0, are similar: 


Akp 45Ti -|- h* 

^ + ssj 7 }' 


(2.15) 


Defining 


m*{N,h*) = 5-f i - 


1 Ak(/i*) 

m ^ Hk(h*) 

k^ko ’ 


2V25V^i(Ji - 2 J 2 ), 


( 2 . 11 ) 


(2.16) 














4 



L 


Figure 2. Regularizing staggered magnetic field h* in the 
J 1 — J 2 antiferromagnet such that SW corrected magnetization 
vanishes. 


we obtain a self-consistent equation for h* 


h* = 4SJi 



(2.17) 


In the limit iV ^ 1, m* —>■ toafj and we have 


2SJi 1 


(2.18) 


As seen below, it is essential to determine the actual value 
of h*{L) with a high precision in order to compute ac¬ 
curately various finite size correlations. Since the field 
h* gets rapidly very small with increasing system sizes, 
we resort to a multiple precision evaluation of the self- 
consistent equation Eq. (2.17). In Fig. 2 we present the 
result showing the behavior of h*{L) for some represen¬ 
tative values of J 2 and S. In all cases, the staggered field 
vanishes very fast and is well described by Eq. (2.18) at 
large enough L. 

Interestingly this small staggered field opens a gap in 
the excitation spectrum 


A* ~ ^/2SJlh* 
_ 2SJi 1 
rriAF N 


(2.19) 


which scales in the same way as the Anderson tower of 
states^^. Therefore, the excitation spectrum has linearly 
dispersing Nambu-Goldstone (SW) modes with a level 
spacing ^ 1/L and a tower of states like finite size gap 
^ 1/L^ produced by the symmetry restoring staggered 
field. 

We use this modified finite-size SW approach to com¬ 
pute the entanglement entropy as detailed below. In or¬ 
der to show that it reproduces fairly well the physics of 
finite-size systems, we also compare in Appendix A re¬ 
sults for the finite-size structure factor for S = 1/2 and 
various J 2 < 0 to the ones obtained with the exact QMC 
method. 


C. Entanglement entropy 


As the diagonalized Hamiltonian Eq. (2.10) is non¬ 
interacting, Wick’s theorem eases the computation of en¬ 
tanglement entropy, which can nicely be extracted from 
the correlation matrix^^, an object which contains all 
two-body correlations within a block of sites. For com¬ 
pleteness, we recapitulate here the essential formulae. 

We first need to define single particle Green’s function 
(blbj) = -^ + h and {b^b^) = gij, with 

_ J_ Akcos [k ■ (rj - rj)] 

•'b' - 27V Hk 

k 

s., = -Tw«if£dM£i:i£d. ( 2 . 20 ) 

2N ^ Hk ^ 

k 

We remark that gij = 0 (/^ = 0) if i and j belong to the 
same (different) sublattice(s). 

The entanglement entropy of a region H containing Nq 
sites can then be extracted^^^^® from the eigenvalues lyf 
of the Nfi X Nq correlation matrix C 

= E(/»' + - 9^'J) (2.21) 


where i,j G H. Due to the sublattice properties of / and 
g, we have that Cij = Cji if i and j belong to the same 
sublattice, Cij = —Cji otherwise. 

The Renyi entanglement entropy is obtained as^^ 


Nn 


91 


which for q = \ reads 
Nn 

5*1 = ^ ^ e -b I) In (^121 + , (2.23) 

1=1 £=±1 

and for q = 00 


Nn 


^00 =51 In 


1=1 


D + 2 I ■ 


(2.24) 


As first shown by Srednicki^®, Gallan and Wilczek^^, the 
entropy of a free massless bosonic field obeys a strict area 
law, which is what we observe (data not shown) in the ab¬ 
sence of the regularizing staggered field h*. However, as 
we will see below, the finite staggered field which opens 
a finite size gap ~ 1/N leads to an additive logarith¬ 
mic correction proportional to the number of Goldstone 
bosons. 


III. RESULTS FOR EE 
A. Strip geometry 

Let us start with the case of an L x strip subsys¬ 
tem embedded in an L x L torus, as depricted in Fig. 1 
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Figure 3. Entanglement Renyi entropies for the strip sub¬ 
system with different aspect ratios tjL (upper panel) and fit 
results for the prefactor Iq of the logarithmic scaling term as 
a function of the minimal system size Lmin included in the fit 
(lower panel). The results displayed here have been obtained 
for S' = I and J 2 = 0. Clearly, Iq = 1 independent of q and 
the aspect ratio of the subsystem. 


plotted in the lower panel of Fig. 3 for various values 
of the Renyi parameter and several aspect ratios. For 
q = 1 , 2 we clearly observe that Iq = 1 over basically the 
whole range of Lmin, whereas for larger values of 9 , the 
convergence is relatively slow as these results are to our 
experience hampered by more severe finite size effects. 
Nevertheless, the resulting Iq is already very close to 1 
and the deviation decreases slowly as Lmin is increased. 
This leads us to the conclusion that, within our SW ap¬ 
proach, we find Iq = ncj^ = 1 to be independent on q 
and the aspect ratio of the subsystem, in perfect agree¬ 
ment with the field theoretical result by MG '. 


B. Square subsystems: corner contributions 


In addition to the breaking of continuous symmetries, 
logarithmic corrections to the area law can also be caused 
by geometry: In particular, logarithmic corrections in¬ 
duced by sharp corners of the subsystem have been 
discussed in several works^^ds. 14,16-18,20-22,25,38,39^ 
prefactor of the logarithmic corner correction term is ex¬ 
pected to be universal for all systems with the same type 
of symmetry breaking/phase transition. However, such 
corrections are quite difficult to capture with QMC since 
the prefactor is very small. Together with the contribu¬ 
tions coming from Nambu-Goldstone modes Eq. (1.4), we 
expect a total correction of the form 


(a). This geometry has no corner and we therefore ex¬ 
pect the expression Eq. (1.3) to hold. Results obtained 
from the exact diagonalization of the correlation ma¬ 
trix C for systems up to ^ 10^ lattice sites are shown 
in the upper panel of Fig. 3 where the Renyi entropies 
for q = 1,2,3 are displayed for three representative as¬ 
pect ratios IjL = 1/2, 1/4, 1/8. Note that for this 
strip geometry, translation symmetry of the subsystem 
is used, allowing the diagonalization procedure to reach 
large sizes. This plot clearly demonstrates the area law 
behavior Sq ^ a^L since the dominant scaling behavior 
does not depend on the number of subsystem sites but 
only on its perimeter 2L, which is independent of the as¬ 
pect ratio of the subsystem. The properties of the area 
law prefactor Uq will be analyzed in detail in Sec. IIIG, 
and the universal additive constant 7 ^ from Eq. (1.3) in 
Sec. IV. 

Here, we want to focus on the logarithmic correction 
associated to the breaking of SU(2) rotational symmetry 
with uq = 2 Nambu-Goldstone modes, expected to be 
^ In L. This correction is believed to be universal as it 
should not depend on the geometry and only reflect the 
nature of the continuous symmetry which is broken in 
the ground state^’^do. Therefore, we perform fits to the 
general scaling ansatz 

Sqi^L) = dqL IqYli L bq -\- CqJ L -\- dq/ L^, (3 - 1 ) 

over various fit ranges [Lmin, Lmax]- Results for Iq are 


( 3 - 2 ) 

f/L = i f/L = i 




Figure 4. Difference of entanglement entropies for the S = 
1/2 J 2 = 0 Heisenberg antiferromagnet of square and strip 
subsystems having the same boundary length. The remaining 
dominant scaling term is the logarithmic term which stems 
from the corners of the square subsystem. We show fits (full 
lines) to the form 8Zg(7r/2) ln(L) + bq + CqjL + dq/L^. SW 
results (symbols) are shown for two different aspect ratios a 
for q — 1, 2,3 and 4. 
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9 = 1 

g = 2 

g = 3 

g = 4 


-0.0118 

-0.0064 

-0.0051 


This work 

-0.0118(1) 

-0.0064(1) 

-0.0050(1) 

-0.0043(1) 


Table I. Prefactor Zg(7r/2) of the corner logarithmic correction 
obtained after fitting data in Fig. 4. A comparison with data 
of Casini and Huerta (CH)^^ is also given. 

where the sum is taken over all sharp corners inside the 
subsystem making an angle tpc- Here we aim at numeri¬ 
cally extracting Z^(7r/2) for a square subsystem (panel (b) 
of Fig. 1), expected to coincide with the result of a free 
scalar field^^. To do so, we work with a. L x L torus and 
substract the entropies of a periodic (corner-free) strip 
of size L X £ from those of a L/2 x L/2 square. Both 
subsystems having the same area law ^ 2L, independent 
of the strip aspect ratio and identical logarithmic 
corrections due to Goldstone modes, we therefore expect 
the leading term of this difference to be given only by the 
corner log contribution: 

S'^qu^e _ 5 'Btrip ^ 8Z^(7r/2) InL -h • • • (3.3) 

Numerical results are plotted in Fig. 4 where we clearly 
see that the above difference Eq. (3.3) is indeed domi¬ 
nated by a logarithmic scaling which allows us to extract 
Small variations of the results for different as¬ 
pects ratios of the strips (see left and right panels of 
Fig. 4) can be used as a measure of the error due to fi¬ 
nite size effects and fitting procedure. Our results are 
displayed in Table I where we compare to the free-held 
results by Casini and Huerta (CH)^^. 


Interestingly, we can also study the dependence on the 
Renyi index for non-integer values of q. In Fig. 5 we 
show 1^{tt(2) versus the Renyi parameter q for four dif¬ 
ferent aspect ratios. For q not too large, the estimates 
obtained after fitting our numerical data (see caption of 
Fig. 4), are clearly independent of the aspect ratio, as 
expected. This non-trivial g-dependence for a free scalar 
field can be compared to recent numerical results for 0(1) 
and 0(2) Wilson-Fisher critical points^^, featuring qual¬ 
itatively similar behaviors. 


C. J2-dependence and area law prefactor 

Besides universal contributions arising from Nambu- 
Goldstone modes and corners, we now study the domi¬ 
nant part which governs the entanglement growth with 
the subsystem area. As already discussed in the begin¬ 
ning of the paper, the Ji — J 2 spin-^ Heisenberg model on 
the square lattice is Neel ordered for J2/J1 < 0.5 in the 
large S limit (see Appendix A and Fig. 15 for the critical 
value of J 2 as a function of S). Scanning across the en¬ 
tire Neel ordered regime, we have performed fits to the 
form Eq. (3.1) for various values of the second neighbor 
coupling J 2 and spin S for the strip geometry (corner- 
free) with a 1/8 aspect ratio. Shown in Fig. 6, the area 
law coefficient aq displays a quite remarkable behavior. 
First, the results appear to be almost independent of the 
spin size S. Then, as expected from the mean-field limit 
JijJ 2 —t —00 (see Refs. 40 and 41 and Appendix B), aq 
goes to zero in the limit —J2/J1 ^ 1. This is because 
the ground-state becomes more and more classical, with 



Figure 5. Logarithmic contribution of a ^ corner of the sub¬ 
system as a function of q. This result is obtained for J 2 = 0 
and 5" = I by subtracting the entanglement entropy of a strip 
subsystem with the same perimeter as the square subsystem 
with 4 ^ corners and fitting to the same form as shown in 
Fig. 4 for different aspect ratios of the strip. Up to slight 
deviations for larger Renyi indices q due to finite size effects, 
the results do not depend on the aspect ratio. 





Figure 6. Top: Area law coefficients a, for g = 1, 2 (extracted 
for an aspect ratio £/L = 1/8) as a function of the next neigh¬ 
bor coupling J 2 /Ji. Approaching the critical point in the frus¬ 
trated regime (J 2 > 0) the area law coefficient grows rapidly. 
Bottom: Prefactor Iq of the log correction in Eq. (1.3) due to 
the two Goldstone modes of the antiferromagnet. The devi¬ 
ation of Iq from 1 close to the critical point in the frustrated 
regime reveals the limitation of the SW approximation. 
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a very low entanglement. However, less is known when 
frustration is turned on, and we observe a rapid growth of 
the area law term when the critical point is approached, a 
feature also observed for the unfrustrated Heisenberg bi- 
layer^"*’. Note that the validity of the SW calculation can 
be questioned when quantum fluctuations become large, 
approaching the critical point J^. Nevertheless, we be¬ 
lieve the results to be under control if 1 — mAF/S ^ 1, a 
condition which can be checked in Fig. 15 of Appendix A. 
Moreover, as long as the logarithmic term in the entan¬ 
glement entropies scaling is still present and equal to one 
(due to the two Nambu-Goldstone modes, as shown in the 
bottom panel of Fig. 6 thus confirming the universality 
across the ordered phase), we believe the SW approxima¬ 
tion correctly captures the behavior of EE. In practice, 
we start to see a deviation of Iq from unity only for the 
very last points at J 2 /J 1 > 0.3 in Fig. 6. 

IV. UNIVERSAL ADDITIVE CONSTANT 7°'^'* 
A. Direct extraction from large-S data 

Following MG', a universal additive constant 
depending on the aspect ratio i/L of a strip, 
appears in the Renyi EE scaling Eq. (1.3). Nevertheless, 
there is also a non-universal term involving the spin stiff¬ 
ness and the SW velocity in Eq. (1.3). It is therefore 
much easier to work in the S ^ 00 limit of the Ji — J 2 
Heisenberg Hamiltonian where Ps and v are known ex¬ 
actly. In such a limit and having shown above that the 
logarithmic prefactor is exactly given by Iq = uq/^ = 1 


(the corner constribution vanishes for this geometry), we 
expect the EE scaling for strips with an aspect ratio IjL 
to be in the S' 1 limit 

Sq = a,L +In +l°q^\i/L). (4.1) 

This large S expression has been used to fit our numer¬ 
ical SW data obtained for S = 100 and J 2 /J 1 = 0,-1. 
Results for the additive constant 7 °''‘^(£/L) are plotted 
in Fig. 7 as a function of the aspect ratio ijL for various 
Renyi parameters q. The agreement with the result ex¬ 
tracted from Ref. 7 is excellent for g = 2. The universal 
character of 7 °''‘^(f/L) is also corroborated by the fact 
that our estimates do not depend on the values of the 
second neighbour coupling J 2 /J 1 = 0,-1; the sole J 2 de¬ 
pendence being given by the two first terms in Eq. (4.1). 

For the S = 1/2 Heisenberg model (J 2 = 0), while we 
have seen above that the logarithmic corrections are per¬ 
fectly captured, a precise determination of the additional 
constant 7 °‘'‘^(^/L) is less obvious. Indeed, as shown in 
Fig. 8 for g = 2, using ps/v = 0.11675 from previous I/S' 
estimates'*^, the y”'* estimates are close but do not agree 
with the S = 100 results. Taking instead the most recent 
QMG estimate*^ for this ratio ps/v = 0.10882(4), the 
agreement is clearly better. We have also checked that re¬ 
sults on 7 °''‘*(£/L) obtained by taking into account higher 
orders in 1/S from Ref. 43 give indeed an improvement 
over the 1/S order, but are not as good as the QMG 
result. 
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Figure 7. Universal additive constant for different 

Renyi indices g as a function of the aspect ratio of the strip 
geometry for the S = 100 Heisenberg Ji — J 2 antiferromagnet 
at J 2 = 0. The inset displays a zoom, showing that our 
result is in perfect agreement with the universal geometric 
constant obtained by MG^ for q = 2. Results for J 2 /J 1 = —1 
at g = 2, shown by red plus signs, agree perfectly with MG 
(black circles) and J 2 = 0 (black x), confirming universality. 


Figure 8 . Results for -S' = | and J 2 = 0 for 72 '^'* obtained from 
fits using rectangular subsystems with fixed aspect ratios. For 
S' = |, we show two sets of results using slightly different 
estimates for Ps/v. Data shown in blue use ps/v = 0.10882(4) 
from the most recent QMG estimate*^, while data shown in 
green use Ps/v = 0.11675 from a 1/S calculation*'*. The use 
of the QMG result leads to a much better agreement with our 
S = 100 data. 














B. Connection to the free scalar field result 

In a (corner free) strip subsystem geometry with a fi¬ 
nite aspect ratio f/L, one expects for gapped free bosons 
with a very large correlation length ^ ^ L the following 
subleading corrections to the area law^: 

A5, = iln(|)+7j^'^'=(^/i), (4.2) 

where 7 ^''®®(f/L) is a universal geometric constant which 
depends non-trivialy on both the Renyi parameter and 
the aspect ratio. By artificially gapping the linear SW 
Hamiltonian Eq. (2.4) with a very small staggered field h, 
the dispersion relation in the vicinity of its two minima 
reads 

fl(k) = y'SSJih + 8S'2Ji(Ji - 2J2)|kP, (4.3) 

thus leading to 

(4.4) 

The correction due to the two minima becomes 

A5, = In + 27 ^®(£/L), (4.5) 

which is used to fit numerical SW results with a very 
small field h = 10 “^® to extract 7 j‘'®®(£/L), shown in 
Fig. 9. 

Quite interestingly, from the above formulation we can 
infer a very simple and direct relation between and 
7 ^''®®. Indeed, in the large S limit the size-dependent 



e/L 


Figure 9. Universal additive constant 7 q'®®(f/L) for free 
bosons plotted against the aspect ratio of the strip subsystem. 
Linear SW results for the 5 = 1/2 Heisenberg antiferromag- 
net at J 2 = 0 are shown together with large-5 estimates of 


staggered field (added to artificially restore zero sublat¬ 
tice magnetization) takes the exact form h*{L) = 
Plugging this into Eq. (4.5) yields 

7”d SU(2) ^27*®®-Lln2, (4.6) 

which agrees with MG^, but only when q = 2. In Fig. 9, 
comparing 7 ^®® to ( 7 ”'^ - ln2)/2 for q = 1,4 gives 

a perfect agreement for a wide range of aspect ratios. 

One can also repeat the same argument for the XY 
model with only one Nambu-Goldstone mode^® to get 

^ordU(l)^^free^5^^2. (4.7) 

The reason for the disagreement between our large S 
approach — expected to become unbiased in the limit 
5 —>■ oo — and the result from Ref. 7 for 7 ”^ is not clear 
at the moment. 


C. Limit of vanishing aspect ratio 

In this section we shed light on the divergent behavior 
of 7 °'''^ for small aspect ratios by calculating 7 °“'^ in the 
extreme limit oi i/L ^ 0 using subsystems with a fixed 
number £ of lines and thus a varying aspect ratio as a 
function of L. In order to achieve this, we work with 
S = 100 at J 2 = 0, and we want to subtract all dominant 
terms, in particular the linear area law contribution aqL. 

Let us therefore start with a study of the dominant 
scaling contribution of Sq by plotting Sq/L vs. 1 /L as 
displayed in Fig. 10. We show the area law behavior 
by plotting S2/L vs. 1 /L and an extrapolation L —>■ 
00 , which guarantees to eliminate all subleading terms. 
The figure shows two sets of curves. In the first one, 
each curve corresponds to subsystems with a constant 
aspect ratio, such that 7 ”^ jg a constant for each curve. 
These curves all yield identical area law prefactors = 
0.190216(1) as expected. 

The second set of curves shows results corresponding 
to a fixed number £ of lines in the subsystem (be. a Gleg 
ladder), which implies that the aspect ratio of the sub¬ 
system is a function of L. The dominant linear prefactor 
ttq found for the Gleg ladder subsystem is different from 
the fixed aspect ratio value, which is approached only for 
£ ^ 1 . The reason for this discrepancy lies in the diver¬ 
gent behavior of when £/L tends to zero and the fact 
that the assumption that the only surviving term in the 
scaling of Sq/L at large sizes is the area law is no longer 
true. In fact, as the aspect ratio of the subsystem changes 
constantly, 7 °"''^ seems to have a contribution that is lin¬ 
ear in the inverse aspect ratio and hence leads to a shift 
or an effectively changed area law prefactor. As a next 
step, we will try to determine this contribution. 

Fig 11 shows our data for as obtained in Sec. IV 
multiplied by the aspect ratio as a function of the as¬ 
pect ratio in order to extract the singular contribution 
72 as the intercept at vanishing aspect ratio. The data 
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Figure 10. S 2 /L vs. 1/L for the strip shaped subsystem with 
different aspect ratios (rainbow curves). SW data obtained for 
5 = 100 and J 2 = 0. The dotted lines show an extrapolation 
for infinite system sizes which demonstrates that for all aspect 
ratios the leading scaling term is indeed an area law. We also 
show results for the line and 2,3,4 and 5 leg ladder subsystems 
(labeled lines) which clearly deviate significantly from this 
behavior. This discrepancy is caused by the singular behavior 
of 7 °“^'^ (and higher subleading terms), as the aspect ratio of 
these subsystems changes constantly with L and runs into the 
divergence of small aspect ratios (see text). 

shows convincing evidence that this contribution indeed 
extrapolates to a nonvanishing value, which we deter¬ 
mine by a cubic fit. With this information, we can now 
decompose in a singular and a regular component: 


Figure 11. Singular contribution 72 of 72 '^‘^(|;) (red) and sin¬ 
gular contribution ri 2 of Cq (black). The dashed line corre¬ 
sponds to cubic fits for the smallest aspect ratios and are 
used to extract estimates of 7 J and , given by the intercept 
at vanishing aspect ratio. 72 corresponds to the contribu¬ 
tion of 72 to the linear scaling of the entanglement entropy of 
fixed width subsystems. 772 is the linear scaling contribution 
stemming from the EE scaling term Cq/L. 

Clearly, for fixed aspect ratio subsystems the terms 7 * 
and rj* become irrelevant for the area law for large system 
sizes. However, for fixed width subsystems, the effective 
linear (in L) scaling coefficient is in fact given by 

al = al+j;/i + rj;/i^. (4.10) 


lim 7^0 = 0. (4.8) 

For completeness, we provide a table of 7 * for other Renyi 
indices q in Tab. II. 

In general, we can assume that other subdominant 
terms show pathologic behavior in the limit of vanish¬ 
ing (and non constant) aspect ratios, i.e. for fixed width 
I subsystems, we will for the moment assume that they 
could produce a total correction to the area law of the 
form r\* ^ in total. The scaling of the EE then reads: 

Sq = {a\+^yi + rfje)L 

+ '^ln(^L) + 7 ”d^ (4.9) 


g 

1 

2 

3 

4 

T'9 

-0.07677 

-0.04579 

-0.03530 

-0.03123 


Table II. Values of 7 , for different Renyi indices. 


We can therefore obtain (in the limit of large L) 7 ™'^ from 
fixed width subsystems by subtracting several terms from 
Sq. Obviously we need to subtract (a^ — 7 */.^) L to elimi¬ 
nate the linear contribution (note how this automatically 
takes care of the unknown terms rj*). 

The second term that we have to subtract from the 
EE is the logarithmic term which is due to the spon¬ 
taneous breaking of SU(2) symmetry. We have argued 
above alongside with several works^^® that its value is 
no/2 = 1 for the case of fixed aspect ratio subsystems 
and shown in Ref. 10 that this is also true for fixed width 
i subsystems, such as the single line with £ = 1 , we there¬ 
fore subtract the term no/2\n{^L), taking also care of 
the constant stemming from Ps/v, that we know with 
great accuracy for the case ol S = 100 at J 2 = 0. 

Remaining subleading terms are expected to die off in 
the limit of £/L —>■ 0 and are therefore unimportant in 
the region of interest. 

In total, for the limit of L —>■ 00 and a fixed width £ of 
the subsystem, we obtain 7 ™'^ through: 
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Figure 12. Data for 72 ''^ as obtained from fits for fixed aspect 
ratio subsystems (red x) shown together with the result of Eq. 
(4.11) for fixed width subsystems. The agreement is excellent, 
even for relatively large aspect ratios which correspond to 
small systems in the fixed width case. The inset shows a 
zoom. Lines are guides to the eye. 


We can now apply Eq. (4.11) to calculate 7 ”'^ in the 
small aspect ratio regime from fixed width subsystems 
of width i. Fig. 12 shows our result in comparison the 
previously obtained values of 7 “*^ from fixed aspect ra¬ 
tio subsystems (strips). SW results for £ = 1 and £ = 2 
are built on an analytical derivation (presented in Ap¬ 
pendix C) obtained exploiting the fully symmetric na¬ 
ture of such subsystems. The perfect agreement of the 
results obtained with different methods and in particular 
the agreement of the results for different £ is a strong 
evidence for the reliability of this result and therefore 
demonstrates also the singular nature of given by 
the singular component 7 *. 


£ 


aq+'y*qle + V*ql£‘^ 

< +7j/£ 

1 

0.156142 

0.155936 

0.144426 

2 

0.170287 

0.170199 

0.167321 

3 

0.176265 

0.176232 

0.174953 

4 

0.179543 

0.179488 

0.178768 

5 

0.181492 

0.181518 

0.181058 


Table III. Comparison of the directly obtained linear scaling 
factor Uq of fixed width subsystems to the result obtained 
using the singular contributions 72 and from subdominant 
terms. The inclusion of rj* significantly improves the result 
and provides numerical evidence for the correctness of Eq. 
(4.10). 

Can higher subleading terms generate corrections to 
the area law coefficient? It is certainly legal to assume 
that pathological behavior in the limit of £/L —>■ 0 is 
not only present in the scaling constant 7 °'''^ but also 
in higher terms, such as CqjL and dq/Lp'. However, for 
them to modify the area law coefficient, they have to di¬ 


verge much faster, i.e. as P' jlP for the case of Cq. In 
order to investigate this possibility, we plot C 2 77 in black 
in Fig. II and observe that a (small) nonzero contribu¬ 
tion to the linear scaling in h is indeed present which we 
call 7^2 (here we will neglect the contributions to r\ from 
even higher terms, which are difficult to access through 
fits to numerical data). Let us finally plug all the infor¬ 
mation together and see if the singular contributions of 
subdominant terms can explain the discrepancy between 
a* and observed in Fig. 10 by comparing in Tab. 
Ill results for as obtained from a direct fit to fixed 
width EEs and for a* + 7 |/£ + 77q/£^- The left column 
shows the total linear scaling prefactor for fixed width 
subsystems as displayed in Fig. 10, while the rightmost 
column shows the fixed aspect ratio linear scaling pref¬ 
actor a* corrected by the singular contribution of 
giving reasonable agreement. The middle column takes 
into account the next subdominant singular contribution 
r]q from the term Cq jL a,s discussed above and reproduces 
the direct fit result to very high accuracy, thus providing 
strong evidence for the correctness of Eq. (4.10). We 
expect that even less dominant terms, such as dq/Lp will 
provide further corrections, which are relevant for small 
widths £ and should account for the remaining discrep¬ 
ancy, these terms are however very small and extremely 
difficult to extract numerically. 


V. DISCUSSIONS AND CONCLUSIONS 


In this work, we have performed a high-precision SW 
study of the Ji — J 2 Heisenberg SU(2) antiferromagnet 
on the square lattice in order to investigate its quan¬ 
tum entanglement properties. Numerical calculations on 
finite size systems have been performed with an artifi¬ 
cial restoration of zero sub-lattice magnetization using a 
small size-dependent staggered field h*{L). Several sit¬ 
uations have been explored, and we have obtained finite 
size scaling results at large enough size such that the var¬ 
ious terms appearing in the entanglement entropies have 
been precisely computed. 

The universal logarithmic correction due to Nambu- 
Goldstone modes associated with the breaking of a con¬ 
tinuous symmetry (SU(2) in the present case) are well 
captured, giving a correction perfectly fitted by ^ InL, 
independent of the Renyi index q. In the case of subsys¬ 
tem having sharp corners, additional (negative) logarith¬ 
mic corrections have been precisely evaluated, in perfect 
agreement with scalar field theory predictions^^. The 
Ji — J 2 model also offers a nice playground where we 
could check universality of the logarithmic term across 
the entire ordered regime but where we could further 
study the non-universal area law part which exhibits a 
non-trivial behavior, with a noticeable growth approach¬ 
ing the critical point in the frustrated regime. In the 
opposite limit of strong ferromagnetic second neighbor 
coupling, the mean-field limit is recovered with a vanish- 
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ing area law term and a smooth crossover to a purely 
logarithmic scaling of the entropies. 

Part of this work was also devoted to the study of 
the additional constant term expected to be univer¬ 
sal for strip subsystems^, only depending on their aspect 
ratio. It then appeared crucial to impose zero sublat¬ 
tice magnetization in the finite size SW theory, yielding 
a unique size-dependent staggered field h*{L) which (i) 
mimics a tower of state gap ~ 1 jL? in the excitation spec¬ 
trum (responsible for the logarithmic correction), and (ii) 
leads to the correct additional geometric constant 
in perfect agreement with MG^, at least for q = 2. K sim¬ 
ple and direct relation with non-interacting bosons was 
also derived. Finally we have precisely investigated the 
limit of vanishing aspect ratios using very large ladder 
subsystems in the limit of finite number of legs, discover¬ 
ing that the geometric constant contains both a regular 
part and a singular component in this limit. Our study is 
concluded by showing that singular components of even 
less dominant terms explain perfectly the discrepancy of 
the area law terms obtained from fixed width vs. fixed 


aspect ratio subsystems. 

Among the potentially interesting future directions, a 
quantitative study of the geometric constant using exact 
Monte Carlo, while very challenging, appears to be a 
very important point in order to test the validity of our 
prediction for q > 2. It may also be interesting to extend 
the present SW approach to other continuous symmetries 
like SU(N) models using modified flavor-wave theory for 
instance. Other geometries or d = 3 are certainly of great 
interest also, with a larger choice of subsystem shapes. 
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Appendix A: Details of spin-wave calculations 

This appendix provides details of spin-wave calcula¬ 
tions which are not crucial for the computation of entan¬ 
glement entropy, but which are nonetheless useful for an 
understanding of the method and its range of validity. 
We also provide a comparison between the finite-size SW 
approach and direct QMC computations for S' = 1/2 for 
the antiferromagnetic structure factor in the ferromag¬ 
netic range of next neighbor coupling J 2 < 0. 


1. Spin-wave spectrum and velocity 



We present in Fig. 13 the spectrum flk = 2-\/A^ — B'^ 
in the direction = ky (obtained from expressions 
Eq. (2.6)) as a function of k^, for different coupling 
strenghts J 2 and for a spin value S = 1/2. The inset rep¬ 
resents the spin-wave velocity, Eq. (2.11), as a function of 
J 2 /J 1 . We see that the velocity vanishes at J 2 /J 1 = 0.5 
where the SW spectrum features a continuous line of min¬ 
ima at fcj; = 0 and ky = 0, as depicted in Fig. 14. 


Figure 14. 2D color map of the SW spectrum at h = 0 for 
J 2 /J 1 = 1/2 and S = 1. A line of minima is visible along the 
fcj; = 0 and ky = 0 directions. 

Ji), above which the SW-corrected order parameter van¬ 
ishes, is also represented in the inset of Fig. 15 as a func¬ 
tion of S where we observe that Tf —f 1/2 when S gets 
larger. 



Figure 13. SW spectrum at h = 0 for various J 2 IJ 1 plotted 
along the = ky direction. Inset: SW velocity Usw. 


Figure 15. SW results for the AF order parameter Eq. (2.12) 
of the Ji — J 2 model on the square lattice for various spin 
sizes S. Inset: Critical frustration (in units of Ji) plotted 
against the spin length S. 


2. Range of non-vanishing staggered magnetization 

a. Antiferromagnetic order parameter 

Eq. (2.12) can be evaluated numerically for different 
values of the spin size S and second neighbor coupling 
strength J 2 /Ji, in order to probe the range of validity 
of the spin-wave approach, which assumes an ordered 
ground-state. This is illustrated in Fig. 15 where the 
AF order parameter is represented, and as expected, 
is clearly enhanced by ferromagnetic diagonal coupling 
J 2 /J 1 < 0 while it decreases towards zero when J 2 /J 1 
approaches 1/2. The critical frustration J| (in units of 


b. Finite size SW: AF structure factor 

To illustrate the interest of using a formulation of 
SW which treats more correctly hnite-size systems, we 
present results for the computation of the staggered 
structure factor per site on finite square lattices L x L = 
N: 

= (Ai) 

ij 

Using Wick’s theorem, all two-spin correlators can be 
computed in terms of the fij and gtj functions dehned in 
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is a good starting point to study ground-state properties 
on finite systems and in particular the finite size scaling 
of the entanglement entropy, as discussed in the main 
text. 


Appendix B: Mean-field limit 


In the limit —J 2 /J 1 ^ 1 one should recover the mean- 
field result obtained for example for the Lieb-Mattis 
model^^’^^. In such a limit, perfect ferromagnetic cor¬ 
relations between spins belonging to the same sublattice 
imply fij = S ioi i ^ j both on the same sublattice 
{gij = 0) and fa = S Antiferromagnetic correla¬ 

tions between opposite sublattices yield 


Figure 16. Staggered structure factor per site s( 7 r, tt) of the 
Ji — J 2 antiferromagnet Eq. (1.5) plotted against the inverse 
system length 1 /L for 3 values of J 2 , with Ji = 1. Symbols 
show T — 0 QMC results, dashed lines are quadratic fits of 
the form Eq. (A3), and the full lines are modified SW results 
using Eq (A2). 


Eq. (2.20) of the main text. Imposing that (S^ • Sj) = 
{SfS^) (because the theory is strictly speaking not rota- 
tionally invariant), we obtain: 

= (4+4) - (A2) 

b 

A quantitative comparison between the above SW ex¬ 
pectation and exact quantum Monte Carlo simulations 
is shown in Fig. 16. Ground-state expectation values for 
s(7r, tt) of the Ji — J 2 Hamiltonian Eq. (1.5) with S' = 1/2 
and J 2 = 0, —1, —5 have been obtained for various 
square lattices LxL using the stochastic series expansion 
algorithm^^. One sees in Fig. 16 that the agreement is 
fairly good, in particular for strong second neighbor fer¬ 
romagnetic coupling J 2 /S 1 = —5. Interestingly, the finite 
size scaling behavior, expected from previous works^®’"*^® 

s(7r,7r) = TOaf + wi/L -I- TO 2 /E^ -I - (A3) 

is very well captured by SW calculations, as visible in 
Table IV where QMC and SW estimates for tuaf, ^ 1 , 
and m 2 are compared and show a good agreement. 


J 2 /J 1 

rriAF 

SW / QMC 

mi 

SW / QMC 

m 2 

SW / QMC 

0 

0.092 / 0.093(1) 

0.55 / 0.60(1) 

0.8 / 0 . 6 ( 1 ) 

-1 

0.175 / 0.167(1) 

0.42 / 0.47(2) 

0.4 / 0.15(9) 

-5 

0.225 / 0.223(1) 

0.25 / 0.26(2) 

0.2 / 0.2(1) 


Table IV. Fit parameters from Eq. (A3). 

The fact that finite size corrections are well captured 
by this modified SW formalism is a confirmation that it 


(Bl) 

ij 

thus leading to gij = S + 1/N for i ^ j on opposite sub¬ 
lattices {fij = 0). Therefore non-zero matrix elements of 
the correlation C (for i and j on the same sublattice) are 
given by 


Q, = ,5(1 -rn) + \ 

= ,5(1 - rp), (B2) 

where tq = Nq/N is the ratio between the number of 
sites inside the sub-system Nq and the total number of 
sites N. The spectrum of the correlation matrix C is 



Nn Nn 


Co 


a 

a 

0, 

Pi 


Figure 17. Mean-field limit for the Renyi entropies Sq. The 
symbols show numerical results for J 2 /J 1 = —10® with differ¬ 
ent geometries and spin lengths (5=1, 10), plotted against 
the number Nn of sites in the subsystem fl. The numerical 
results (symbols) are compared to the analytical expression 
Eq. (B5) with rn = 3/8 (a) or rn = 1/2 (b), shown by the 
full lines. 




























































then straigthforwardly given by 


Appendix C: Analytical derivation for 
one-dimensional subsystems 
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^ 1 % = ^5(1 - rn) + ^ (B3) 

^l = \, {l = 3,...,Nn). (B4) 

one sees that only two eigenvalues contribute, in a macro¬ 
scopic way. We then compute directly the Renyi en¬ 
tropies for any partition Nq and any q > 1 


Sg = \n{Nn) + In 


S{1 - rg)' 
2 


+ 2 


Inq 


(B5) 


The area law term vanishes, and the dominant scaling is 
now a pure logarithm of the number of sites JVq. This 
exact expression can be compared to the numerical so¬ 
lution of the SW Hamiltonian for very large negative 
values of J 2 . In Fig. 17 we show numerical results for 
J 2 /Ji = —100000 for two values of S and two different 
geometries, which compares extremely well with the MF 
limit expression Eq. (B5). Note that the lines are not 
fitting functions. If we try instead to fit to the general 
form ttqL + IqhiL -\- hq + CqjL + dq/L'^, we end up with 
Gq ~ 10“^^ and Iq = 2. 


A great simplification for the computation of entangle¬ 
ment entropy is possible if all sites i and j inside a subsys¬ 
tem are equivalent, or in other words if the matrix ele¬ 
ments Cij only depend on the relative distance jr^ — rj|. 
In such a case, for sites on different sublattices Cij = 0. 
This situation is achieved for one-dimensional subsystems 
with one or two lines (Fig. 1 (a) with £ = 1, 2). In these 
specific situations, we can derive analytic expressions for 
the eigenvalues of C, avoiding a numerical diagonaliza- 
tion. This has first been discussed in Ref. 10, and we 
provide here details of this calculation, starting with the 
case of a line-shaped subsystem. 

This subsystem being invariant under translations 
along the x direction, the functions fij and gij defined in 
Eqs.(2.20) only depend on the distance x = Xt — Xj along 
the subsystem. They reduce consequently to 

kx 

^ X! (^os{k^x)(3k, (Cl) 



Figure 18. Area law coefficient aq as a function of the normal¬ 
ized order parameter S — rriAF, which is small in the ordered 
phase. Here, we show results for 5" = The inset depicts 
the prefactor of the logarithmic entanglement entropy scaling 
term, which has a plateau at = 1 for intermediate J2 and 
evolves to the mean field limit of J2 — >■ — 00 . It is unclear if 
the behavior in the crossover region is a finite size effect or an 
artefact of the spin wave method. 

Before concluding, we want to briefly comment on the 
crossover to the MF limit when the ferromagnetic second 
neighbour is turned on towards very large values. This 
is illustrated in Fig. 18 where the rapid decrease of the 
area law coefficient aq (for q = 1,2) is shown versus the 
quantum depletion of the AF order parameter S — toaf- 
In the same time, the log coefficients h and I 2 , plotted in 
the inset of Fig. 18, crossover from Zg = 1 up to Z, = 2 in 
the limit of vanishing quantum fluctuations toaf —f S. 


with 

and ^ ^ (C2) 

Ky Kty 

which satisfy the property 

Since the functions fx and gx possess translation and 
reflection symmetries 

fx — Jl—x — Jl+x: and gx — gh—x — gL+xt (C^) 

so does the correlation matrix: Cij = C{1 = \xi — Xj\) = 
C{L — 1) = C{1 — L). Since furthermore C{x) vanishes 
for odd distances, it is convenient to re-index all sites on 
one sublattice from 1 to L/2 (say, blue sites in Fig. la for 
£ = 1), and sites on the other sublattices from L/2-I- 1 to 
L (say, orange sites in Fig. 1 a) to block-diagonalize C 
onto two identical blocks of size L/2 x L/2. The transla¬ 
tion invariance ensures that each block is circulant, with 
matrix elements C{1) = Ex even“ J2x odd 9x9x-i- 
The eigenvalues lyf of C are given by the properties of 
circulant matrices'*^^: 

/T\ r^/4—1] \ 

vl = c(0) -b (-l)'c j + ^ c(2f) cos ’ 

Z e {0, 1 ,... ^ - 1}, c = 0 if ^mod2 = 1. 

(C4) 

We can even simplify calculations by noticing that fx 
and gx are discrete Fourier transforms of and re¬ 
spectively. Using the convolution theorem on C'(Z), we 
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arrive at the final expression for the L eigenvalues of C can now introduce the discrete Fourier tranforms of the 
of the single-line subsystem: newcomers f{x, 1) and g(x, 1) 


d = (C5) 

with q e {—TT + . tt}. 

A very similar reasoning can be applied to the case 
of a 2-line (ladder) subsystem with 2L sites (£ = 2 in 
Fig. la). It is convenient to re-index sites by labelling 
(in a zig-zag fashion) all (say, blue in in Fig. la) sites 
of one sublattice from 1 to L, and (orange in Fig. la) 
sites from the other sublattice from L -|- 1 to 2L. Again, 
C is block-diagonal with identical circulants blocks with 
matrix elements 0(1) = fx ft-i “ Ex atdt-i with 

fx = /(a;,0) -f f{x,l) and g+ = g(x,0) -h g(x,l). We 


a/c 


-E 


Ak cos{ky) 

rik 


~ B]g_cos{ky) 

and 


(C6) 

to again be able to apply the convolution theorem. We 
finally obtain that the L eigenvalues of one block of C for 
the ladder subsystem are given by: 


2 _ 

^9 “ 4/V 


{ay - ayf - iPg - ^ yf 


(C7) 


with q G {— TT -I- ^,...7r}. Since C has two identical 
blocks, the 2L eigenvalues for the ladder subsystem are 
obtained by doubling the above spectrum. 






